Self-optimizing analysis window sizing method

ABSTRACT

A method of selecting an optimal window length in a digital processing system includes receiving a digital signal and analyzing the signal with a group of Hanning windows having different sizes. The windows being arranged so that they can be scaled to be comparable. The digital signal is windowed with each selected window and a transform is computed. The transformed signal is scaled and corrected. A metric is computed from the resulting signal for each window. The metrics are compared and the window size is selected based on agreement with a user defined metric. The specific window function, shift and scaling are such that the resulting analysis is mathematically equivalent across different window sizes.

STATEMENT OF GOVERNMENT INTEREST

The invention described herein may be manufactured and used by or forthe Government of the United States of America for governmental purposeswithout the payment of any royalties thereon or therefor.

CROSS REFERENCE TO OTHER PATENT APPLICATIONS

None.

BACKGROUND OF THE INVENTION

(1) Field of the Invention

The current invention relates to a method of signal analysis for signalsof unknown duration or bandwidth.

(2) Description of the Prior Art

The formation of overlapped Hanning-weighted analysis windows is a verycommon first stage of processing in many time-series analysisapplications. The FFT typically follows the window analysis. This is thecase, for example, in the short-time Fourier transform (STFT) which is awidely used front-end for a variety of applications including automaticspeech recognition (ASR). It is often necessary to arbitrarily choose ananalysis window size that is a compromise between the desire for goodtime-domain and frequency-domain resolution. In order to handle a widerange of input time-scales, it is sometimes necessary to use multipleanalysis window sizes in parallel. But this approach creates a newproblem—the resolution of ambiguous results. In other words, one has todecide which FFT size is appropriate for a given input data record. Arelated technique is spectrogram combining where STFTs are combined insuch a way that the choice of best FFT size is made at each grid pointin the time and frequency plane. An important step toward solving thecomparison problem is the assurance that the various SIFTrepresentations are in some way “equivalent”. One possible definition of“equivalent” is the existence of an orthonormal linear transformationrelating one analysis at a given window size to another at a differentwindow size. This definition has relevance from a number of differentviewpoints—from linear subspace analysis to statistical methodsincluding the PDF projection theorem. See “The PDF Projection Theoremand the Class-Specific Method”, IEEE Transactions on Signal Processing,Vol. 51, No. 3 (March 2003). There is a trivial case where two analysesare related by a permutation and thus equivalent. Consider a time-serieswhere the total number of samples T is divisible by N₁ and N₂. If weanalyze with a rectangular window function and use no overlap betweenprocessing windows, the two analyses with window sizes N₁ and N₂ arepermutations of the same input samples. Unfortunately, when anon-rectangular window function is used and the processing windows areoverlapped, there is in general no permutation or orthonormaltransformation relating the two analyses.

Thus, there is a need for a method that will allow comparisons betweenwindowed data having different lengths when non-rectangular windowfunctions are used.

SUMMARY OF THE INVENTION

A first object of the invention is to determine the most appropriateanalysis window sizes for an incoming signal.

A second object is to provide a method of analysis for a signal havingpulses of unknown duration.

Accordingly, there is provided a method of selecting an optimal windowlength in a digital processing system includes receiving a digitalsignal and analyzing the signal with a group of Hanning windows havingdifferent sizes. The windows being arranged so that they can be scaledto be comparable. The digital signal is windowed with each selectedwindow and a transform is computed. The transformed signal is scaled andcorrected. A metric is computed from the resulting signal for eachwindow. The metrics are compared and the window size is selected basedon agreement with a user defined metric. The specific window function,shift and scaling are such that the resulting analysis is mathematicallyequivalent across different window sizes.

Other objects and advantages of the present invention will becomeapparent from the following description, drawings and claims.

BRIEF DESCRIPTION OF THE DRAWINGS

A more complete understanding of the invention and many of the attendantadvantages thereto will be readily appreciated as the same becomesbetter understood by reference to the following detailed descriptionwhen considered in conjunction with the accompanying drawings whereinlike reference numerals and symbols designate identical or correspondingparts throughout the several views and wherein:

FIG. 1 is a diagram showing apparatus that could be used to practice theembodiments; and

FIG. 2 is a block diagram showing one embodiment of the invention.

DETAILED DESCRIPTION OF THE INVENTION

FIG. 1 provides an apparatus 10 for executing this method. An unknownenvironmental signal 12 is received at a sensor 14. Analog data fromsensor 14 is converted into digital data by analog to digital convertor16 at a sample rate. Digital data from analog to digital convertor 16 isjoined to a signal processor 18. Signal processor 18 is capable ofapplying at least one windowing function 20 to the digital data forselecting a time or frequency resolution for analysis. If multiplewindowing functions are applied, the signal processor can apply thesewindowing functions sequentially or concurrently. The windowed data isprovided to a linear transform 22 such as a fast Fourier transform,discrete cosine transform, discrete wavelet transform or the like toproduce a transform. The transform provides the amplitude of the signalaccording to the various basis functions of the transform. In the caseof the FFT, this means according to various frequencies. These groups ofamplitudes are known as bins. Parameters related to the bins can beanalyzed in an analysis component 24 to give the characteristics of theenvironmental signal as user output 26 or as input to another process.

The method of this invention is particularly designed for acoustic data;however it can also be applied to other types of data, such as thatassociated with seismic, radar and radio signals.

The mathematics behind this method is given in the following text. Inputdata sample X has a length T such that X=[x₁, x₂ . . . x_(T)]′. Length Tis given in a number of digital samples. It is desirable to analyze Xinto length N overlapped, shaded windows. The windows are shifted by Ksamples at each update. D is the ratio of window size to time shift suchthat:

$\begin{matrix}{D = \frac{N}{K}} & (1)\end{matrix}$

With T being a multiple of K and X being circularly indexed so thatx_(T+i)=x_(i), there will be exactly DT/N processing windows. The lengthN processing windows are denoted by:y _(t) =[x _(K(t−1)+1) w ₁ ,x _(K(t−1)+2) w ₂ , . . . x _((K(t−1)+N) w_(N)]′  (2)where N is the window length and w={w₁, w₂ . . . w_(N)} is the windowfunction. Y is the complete window analysis, the concatenation of thevarious windowed analysis windows, Y={y₁, y₂ . . . y_(DT/N)}. Theconversion from X to Y is a linear transformation that can be written asY=A X where A is a DT×T matrix and Y is the length DT concatenation ofthe DT/N window y_(t).

Two critical properties occur when using the Hanning weight function.The Hanning weight function is given as:

$\begin{matrix}{{w_{i} = \frac{1 + {\cos\left( \frac{2\;{\pi\left( {i - 1} \right)}}{N} \right)}}{c}},{1 \leq i \leq N}} & (3)\end{matrix}$where c is a constant.

The overlap—add reconstruction of the input X is possible if the shiftedwindow functions add to a constant. This means that the rows of matrix Amust add to a constant. The sum of row i of A is:

$\begin{matrix}{{{\sum\limits_{k = 0}^{D - 1}\;\frac{w_{i + k}N}{D}} = {\frac{D + {\sum\limits_{k = 0}^{D - 1}\;{\cos\left( \frac{2k\;\pi}{D} \right)}}}{c} = \frac{D}{c}}};{D \geq 2}} & (4)\end{matrix}$wherein the summation in the numerator resolves to 0 when D≧2. Anotherrequired property is that the square of the shifted window functionsmust add to 1 or that the squares of the elements in each row of matrixA must add to 1. This does not happen for Hanning-2 (D=2) weightings,but it does happen for all weightings with D≧3 where the sum of thesquares of row i of A is:

$\begin{matrix}{{{\sum\limits_{k = 0}^{D - 1}\; w_{i + \frac{kN}{D}}^{2}} = {{\sum\limits_{k = 0}^{D - 1}\;\left\lbrack \frac{1 + {\cos\left( \frac{2k\;\pi}{D} \right)}}{c} \right\rbrack^{2}} = {{\frac{D}{c^{2}} + {\sum\limits_{k = 0}^{D - 1}\;\frac{\cos^{2}\left( \frac{2{k\pi}}{D} \right)}{c}}} = {{\frac{D}{c^{2}} + \frac{D}{2\; c^{2}}}\; = \frac{3D}{2\; c^{2}}}}}}{{if}:}} & (5) \\{c = \sqrt{\frac{3D}{2}}} & (6)\end{matrix}$the result is 1. It follows, as it explained below, that there is alinear matrix that will allow comparison between different window sizes.

By these principles, there is a full-rank orthonormal transformationQ^(N) ₁ ^(N) ₂ that allows transformation between the window analysis atN₁ and the window analysis at N₂. Mathematically this is given as:Y ^(N) ² =Q ^(N) ¹ ^(N) ² Y ^(N) ¹   (7)

As described above, N₁ and N₂ must both be divisible by D and D≧3. Thenumber of samples T must also be divisible by N₁/D and N₂/D.

Under the conditions given above, the columns of A are orthonormal,therefore:A′A=I  (8)where I is the T×T identity matrix. Equation (8) is a direct result ofthe fact that the row sums of A are constant and the sum of the squaredrow elements equal 1. Using the matrix form of the overlap-add method, Xcan be restored from Y^(N) ₁ using:X=(A ^(N) ¹ )′Y ^(N) ¹   (9)

Y^(N) ₂ can be generated using:Y ^(N) ² =A ^(N) ² X=A ^(N) ² (A ^(N) ¹ )′Y ^(N) ¹   (10)The matrix{tilde over (Q)} ^(N) ¹ ^(N) ² =A ^(N) ² (A ^(N) ¹ )′  (11)is a DT×DT linear transformation of rank T that converts from windowanalysis Y^(N) ₁ to Y^(N) ₂, so it is rank-deficient. Note that A^(N) isrank T. The missing rank can be restored using the orthogonal subspaceof rank (D−1)T. The extended DT×DT matrix is:Q ^(N) =[A ^(N) Ã ^(N)]′  (12)where Ã^(N) is any DT×(D−1)T matrix of orthonormal columns spanning thesubspace orthogonal to the columns of A^(N). Since:Q ^(N) ¹ ^(N) ² =Q ^(N) ¹ Q ^(N) ² ′=A ^(N) ² (A ^(N) ¹ )′+Ã ^(N) ² (Ã^(N) ¹ )′  (13)then it follows that(Q ^(N) ¹ ^(N) ² )′Q ^(N) ¹ ^(N) ² =I.  (14)

The Hanning-D analyses can be viewed purely in terms of orthonormalmatrices or rotations. A linear orthogonal transform such as the fastFourier transform (FFT), discrete cosine transform (DCT) or discretewavelet transform (DWT) can occur after this rotation. This transformisn't required to be the same in each branch and can also be implementedas an orthonormal matrix multiplication. Thus, all points can betransformed by a series of coordinate system rotations.

FIG. 2 provides a block diagram showing application of this method. Adigital signal is provided to the process in block 30. Digital signal 30is provided to a plurality of functions 32 as described above. These canbe a set of Hanning-3 functions having sample lengths of N₁, N₂, . . .N_(I) where I is the total number of functions 32. The set of samplelengths can be calculate geometrically, linearly or at values determinedby a user based on an estimate of the components of the input signal.The output of each of these functions 32 is provided to a transform 34of the same length, N_(i). For most purposes a fast Fourier transformcould be used, but this transform 34 could be other linear transformssuch as a discrete cosine transform, discrete wavelet transform or thelike.

The transformed output is scaled in block 36 by the factor

$\sqrt{\frac{D}{{Ni}^{\prime}}}$the square root of the ratio of D to the length N_(i) of the transform.

A statistic is computed for each length N_(i) in block 40. Thisstatistic summarizes the useful information in the transform outputbins. An example of a statistic is to compute the mean square binamplitude, and/or to locate the largest bin amplitude. Becauseindependent of the window size N_(i), there are the same number DT oftotal bins, and because of the scaling in block 36, the bins have thesame statistics under noise only input. Therefore, it is possible tocreate a statistic in block 40 that has the same output statistics undernoise only hypothesis no matter what the window size N_(i). As providedabove, this statistic is comparable among the different windowingfunction lengths and transforms because of the windowing functionparameters and the scaling. A preferred length N_(i) can then beselected from this statistic in block 42. This preferred length can beutilized either to give information about the signal such as pulselength or to select the best linear transform type and transform lengthfor the given signal.

A first example of using this method is in detecting pulses of unknownduration. Searching for the optimal Fast Fourier Transform (FFT) size ona given digitized input signal is tantamount to pulse length estimation.This method was demonstrated utilizing synthetic data that consisted of768 samples of independent, identically distributed, Gaussian noise plusa short Hanning-weighted sinusoidal pulse of random frequency andstarting time. The pulse length of the short pulse, denoted by P, was96, 192, or 384 samples. The signal to noise ratio was quite low, andthe pulse was difficult to find visually in the time-series. In order tomaintain detectability at short pulse lengths, the amplitude was variedto the P-dependent value 15P^(−1/2) at the peak of the Hanning-shapedenvelope.

The following procedure was performed to form a detection statistic thathas identical noise-only distribution independent from the analysiswindow size. A set of sizes was determined as being likely window sizes.These include the following number of samples, N: (72, 96, 144, 192,288, 384, 576, 768). For each of these sizes, a Hanning-3 analysis (D=3)was performed. The fast Fourier transform of the analysis window iscalculated for each number of samples. The resulting bins are scaled by

$\sqrt{\frac{D}{{Ni}^{\prime}}}$giving a mean square bin value of 1 under the statistical hypothesisthat the input is independent Gaussian noise of variance 1 (called theH₀ hypothesis). Next, the magnitude-square of each bin is calculated.Because the transform is the FFT in which the two end bins arereal-valued, whereas the rest are complex valued, we need to speciallyhandle the end bins. Therefore, the two end magnitude-squared bins(representing the zero frequency and Nyquist frequency) are addedtogether and divided by two, giving an end-bin value. The end-bin valueand the remaining N/2−1 magnitude-squared bins represent a set of N/2exponentially distributed bins having mean 1 under H₀, satisfying therequirement that the statistics are independent of the transform size Nunder H₀.

Once the exponentially-distributed bins are computed for all theanalysis windows for the chosen analysis window size N, theexponentially distributed bins of all the analysis windows are combinedinto a vector w which will have a constant length of DT/2 regardless ofN. In this example, this results in (3×768/2)=1152 exponentiallydistributed bins regardless of N. As a summarizing statistic, the powerlaw statistic is computed for each analysis window size according to:

$\begin{matrix}{z_{N} = {\sum\limits_{i = 1}^{{DT}/2}\;{w\;}_{i}^{\gamma}}} & (15)\end{matrix}$where N is the chosen analysis window size and γ=2.5. This is known asNuttall's power law detector.

Under the H₀ hypothesis, the distribution of w does not dependsignificantly on N. It is always a set of DT/2 bins having standardexponential distribution with mean 1. There is, however, statisticaldependence between bins that is imparted on w by the Hanning weighting,but this is not important if the elements of w are treated as unordered.Under the H₀ hypothesis, the distribution of the feature z_(N) must beessentially independent of the number of samples, N.

In the demonstration, two hundred independent time-series were producedat each pulse length. Each sample was analyzed with each analysis windowsize N and z_(N) was evaluated. The pulse length estimate was:{circumflex over (P)}=arg max_(N) z _(N)  (16)In experiments, it can be seen that most of the estimates are at or nearthe true pulse length, while for noise only, they are widelydistributed.

The foregoing description of the preferred embodiments of the inventionhas been presented for purposes of illustration and description only. Itis not intended to be exhaustive nor to limit the invention to theprecise form disclosed; and obviously many modifications and variationsare possible in light of the above teaching. Such modifications andvariations that may be apparent to a person skilled in the art areintended to be included within the scope of this invention as defined bythe accompanying claims.

What is claimed is:
 1. A method of selecting a transform length for useon a digital processing system comprising the steps of: receivingdigital data related to a periodic physical signal as a plurality ofsamples taken over time; selecting a maximum, linear transform having anoverlap factor D that is at least three and a size N_(I) that isdivisible by D; selecting a summarizing statistic based on user criteriaconcerning transform length such that the summarizing statistic isindependent of window size N_(i) under a noise-only hypothesis;determining a family of Hanning windowing functions of lengths N₀, N₁, .. . N_(I), the length of each windowing function being identified asN_(i), each windowing function shifted by K=N_(i)/D samples andoverlapped by N_(i)−N_(i)/D samples; conducting a set of analyses of theinput signal utilizing each member N_(i) of the determined family ofHanning windowing functions, each analysis providing a windowed signalat the window length N_(i); calculating the linear transform inaccordance with the selected linear transform of length N_(i) eachwindowed signal N_(i), said calculated transform having a plurality ofbins, each bin representing the amplitude of the windowed signalaccording to a particular basis function of the selected lineartransform; scaling the bins from the calculated transform of eachwindowed signal to make the bins independent of window size N_(i);computing the summarizing statistic from the scaled bins of eachwindowing function window size N_(i) in a way that is independent of thewindow size N_(i) under a noise-only hypothesis; comparing saidsummarizing statistics of each windowing function window size N_(i) todetermine the transform size N_(i) that is closest to said user criteriaconcerning transfer length; and utilizing the determined transform sizeN_(i) in analysis of said received digital data.
 2. The method of claim1 wherein lengths N₀, N₁, . . . N_(I), of said family of Hanningwindowing functions are geometrically distributed throughout the rangeof lengths N₀ to N_(I).
 3. The method of claim 1 wherein lengths N₀, N₁,. . . N_(I), of said family of Hanning windowing functions are linearlydistributed throughout the range of lengths N₀ to N_(I).
 4. The methodof claim 1 further comprising steps of: estimating potential signalcomponents from the input signal; and determining lengths N₀, N₁, . . .N_(I), of said family of Hanning windowing functions based on saidestimated potential signal components.
 5. The method of claim 1 whereinthe step of scaling comprises scaling each bin by the square root of theproduct of the number D and the reciprocal of the window length, N_(i).6. The method of claim 1 wherein: said selected linear transform is afast Fourier transform; said summarizing statistic is a power lawdetector giving a power law statistic value for each windowing functionwindow size N_(i); and said step of comparing includes determining thetransform size N_(i) that has the maximum power law statistic value andselecting that transform size N_(i) as the transform size being closedto said user criteria.
 7. The method of claim 1 wherein said step ofutilizing the determined transform size comprises utilizing thedetermined transform size as an estimate of pulse lengths.
 8. The methodof claim 1 further comprising the step of transforming said receiveddigital data into the frequency domain utilizing the determinedtransform size N_(i).
 9. The method of claim 1 wherein said periodicphysical signal is an audio signal.
 10. The method of claim 1 whereinsaid periodic physical signal is a radio signal.
 11. A method ofselecting a transform length for use on a digital processing systemcomprising the steps of: receiving digital data related to a periodicphysical signal as a plurality of samples taken over time; selecting anoverlap factor D as a whole number of at least three; selecting a familyof window sizes N₀, N₁, . . . N_(I), the length of each window sizebeing identified as N_(i), such that each N_(i) is equal to a time shiftK multiplied by the overlap factor D; selecting linear transforms oflengths N₁, . . . N_(i), . . . N_(I); selecting a summarizing statisticbased on user criteria concerning transform length such that thesummarizing statistic is independent of window size N_(i) under anoise-only hypothesis; determining a family of Hanning windowingfunctions of lengths N₀, N₁, . . . N_(I), each windowing functionshifted by K=N_(i)/D samples and overlapped by N_(i)−N_(i)/D samples;conducting a set of analyses of the input signal utilizing each memberN_(i) of the determined family of Hanning windowing functions, eachanalysis providing a windowed signal at the window length N_(i);calculating the linear transform in accordance with the selected lineartransform of length N_(i) each windowed signal N_(i), said calculatedtransform having a plurality of bins, each bin representing theamplitude of the windowed signal according to a particular basisfunction of the selected linear transform; scaling the bins from thecalculated transform of each windowed signal to make the binsindependent of window size N_(i); computing the summarizing statisticfrom the scaled bins of each windowing function window size N_(i) in away that is independent of the window size N_(i) under a noise-onlyhypothesis; comparing said summarizing statistics of each windowingfunction window size N_(i) to determine the transform size N_(i) that isclosest to said user criteria concerning transfer length; and utilizingthe determined transform size N_(i) in analysis of said received digitaldata.
 12. The method of claim 11 wherein said step of utilizing thedetermined transform size comprises utilizing the determined transformsize as an estimate of pulse lengths.
 13. The method of claim 11 furthercomprising the step of transforming said received digital data into thefrequency domain utilizing the determined transform size N_(i).
 14. Themethod of claim 11 wherein said periodic physical signal is an audiosignal.
 15. The method of claim 11 wherein said periodic physical signalis a radio signal.